Shape changes and motion of a vesicle in a fluid using a lattice Boltzmann model 
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We study the deformation and motion of an erythrocyte in fluid flows via a lattice Boltzmann method. To 
this purpose, the bending rigidity and the elastic modulus of isotropic dilation are introduced and incorporated 
with the lattice Boltzmann simulation, and the membrane-flow interactions on both sides of the membrane are 
carefully examined. We find that the static biconcave shape of an erythrocyte is quite stable and can effectively 
resist the pathological changes on their membrane. Further, our simulation results show that in shear flow, the 
erythrocyte will be highly flattened and undergo tank tread-like motion. This phenomenon has been observed 
by experiment very long time ago, but it has feazed the boundary integral and singularity methods up to the 
present. Because of its intrinsically parallel dynamics, this lattice Boltzmann method is expected to find wide 
applications for both single and multi-vesicles suspension as well as complex open membranes in various fluid 
flows for a wide range of Reynolds numbers. 
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Vesicle whose membrane consists of lipid bilayer is es- 
sential to the function of biological systems ll]]. Erythro- 
cyte is the most important kind of vesicles. In the past 30 
years, the dynamics of vesicle has received particular atten- 
tion 121 |3l |4l |5[]. It has been recognized that the equiUbrium 
shape can be obtained by minimizing the bending energy of 
the membrane iQl- However, the studies on the unsteady states 
lag behind, mainly because of the numerical difficulty on cap- 
turing the coupling between the membrane and ambient fluids 
while the vesicle is deforming and moving under the hydro- 
dynamic forces exerted on both sides of its elastic membrane 
yD. When the vesicle is very close to a solid static bound- 
ary and the Reynolds number is very small, lubrication theory 
has been extended to the study 101 ■ Recently, the deforma- 
tions of vesicles in the approximation of Stokes flow has been 
extensively studied by the boundary integral and singularity 
methods (|3l|5|], which has a limitation imputed to the absence 
of numerical smoothing and accuracy of regridding so that 
it can not describe the highly flattened tank-treading shapes. 

Erythrocytes in large blood vessels, which have larger 
Reynolds numbers, not only affect the viscosity of the 
fluid 101, but also often subject to pathological changes of 
their membranes due to the large shear stress [80. This calls 
for computing models for deformations of vesicles, partic- 
ularly with inhomogeneous membrane, in fluid flow with a 
wide range of Reynolds numbers. Moreover, considering the 
numerical complexity of the vesicle deformations and multi- 
vesicle suspensions, and the recent development of compu- 
tational technique, especially on the PC clusters and internet 
grid computing, it is most desirable that the numerical mod- 
els are intrinsically parallel. The lattice Boltzmann method 



(LBM) Jgl [l^l has localized kinetic nature which is not only 
intrinsically parallel but also easy to capture the interaction 
between a fluid and a small segment of a deformable bound- 
ary. In the past fifteen years, the LBM has been recog- 
nized as an alternate method for computational fluid dynam- 
ics, especially in the areas of complex fluids such as par- 
ticle suspension flow ll ill , binary mixture L121 and blood 
flow ill3ll . However, applicability of LBM on the vesicle de- 
forming and moving remains an open problem. A crucial diffi- 
culty is the physics that the membrane is viscoelstic and area- 
unchangeable while the total volume of the vesicle keeps a 
constant. Numerically, the viscoelstic property demands that 
the boundary condition and the hydrodynamic force should 
be accurately treated even in a very small segment. In two- 
dimensional problems, area- and volume- unchangeability are 
reduced to inextensibility of the one-dimensional membrane 
and area-unchangablility of the area surrounded by the one- 
dimensional membrane. In this Letter we develop a lattice 
Boltzmann model in which the elastic properties, length con- 
strain of the membrane and volume constrain on area of vesi- 
cle incorporated in a discrete way with a boundary treatment 
for moving boundaries. 

We choose to work with the D2Q9 model on a two- 
dimensional square lattice with nine velocities ifioll . Let 
fa (x, t) be the non-negative distribution function which can 
be thought as the number of fluid particles at site x, time t, 
and possessing one of the nine velocities Gq,. The distribu- 
tion function evolves according to a Boltzmann equation that 
is discrete in both space and time, 

/„(x + e„,f + 1) - /«(x,i) - --{fa - f^). (1) 

T 

The density p and macroscopic velocity u are defined by 
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Here, the equilibrium distribution function /^"^ depends only 
on the local density p and flow velocity u. The pressure and 
the viscosity are defined by the equations p ^ c^p with = 
1/3 and = (2t - l)/6, respectively IToll . 

The membrane of the erythrocytes has bending rigidity po- 
tential, which can be written as ||2,|5t] 
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k^dl, 



(3) 



where ks is the bending modulus, k and I are the curvature 
and the arc length of the membrane, separately. Biomem- 
branes are formed by a lipid bilayer, which is viscoelastic. The 
erythrocyte viscoelasticity is usually assumed to be Kelvin- 
Voigt 1 141] and described by 



T 



(4) 



where T,„„ is the viscous stress, and ry are the strain rate 
and the viscous coefficient of the membrane. The viscoelas- 
ticity of the membrane does not change the static shapes of 
erythrocytes. 

Numerically, the membrane of a two-dimensional erythro- 
cyte is discretized into equilength segments. We implement 
a no-slip boundary condition and compute the hydrodynamic 
forces on both sides of each segment according to the scheme 
proposed recently llSll . The isotropic dilation can be de- 
scribed by adopting an in-plane potential 0^ between neigh- 
boring segments as 



1 ^ 



(5) 



where is the elastic coefficient of the membrane, Iq and 1^ 
are the original and simulated length of segment i respectively. 
Pozrikidis has already shown that the transverse shear tension 
due to the bending energy can be computed as Isfl 



F^ = kH^ 



dk 



(6) 



The force exerted on segment i on the normal direction due to 
the bending energy is 



(7) 



where /c^+i and ki are the curvatures of the membrane at seg- 
ments I + 1 and i, respectively. The membrane viscous resis- 
tance on the normal direction on segment i is 



(8) 



where and u,;,„ is the velocity of segments i + \ and i 

along the normal direction of segment i, separately. 

The translation of each segment is updated at each New- 
tonian dynamics time step according to the sum of all the 
forces on the seg ment by using a so-called half-step 'leap- 
frog' scheme lll6ll . 



The membrane parameters are set to be fcs = 1.8 x 10 
dyn • cm Jstl and -q = 1.0 x 10~^ dyn • s/cm 1 14|. The blood 



serum is usually assumed to be Newtonian and has a viscosity 
V ~ 0.01 cm^/s and density p = 1.00 g/cvn^ fs]]. The fluid 
inside the erythrocytes is also assumed to be serum too. The 
thickness and the density of the membrane are set to be 0.02 
/im and 1.00 g/cm"^, respectively ITtIi . The cross-membrane 
pressure drop of an erythrocyte can be expressed by the chem- 
ical potential drop jlSll 

/Pout \ 



= i?rin(- 



(9) 



where pout and pin are the pressure outside and inside the ery- 
throcyte, separately. The temperature is set to the human body 
temperature 37° C. The longer axis of an erythrocyte is set to 
3.9 (Um and the according shorter axis is 0.4 p m UtIi . The 
elastic modulus of isotropic dilation is fcfc — 500 dyn/cm Isj] 



FIG. 1: The static profiles of an erythrocyte for A/i = 2.71 (□), 
2.97 (O), 3.39 ( ), and 3.82 (A) x 10"^ J/mol calculated from lattice 
Boltzmann simulations (symbols) together with those from a shoot- 
ing method IT^ (lines), x and y are normalized by the total length 
of the membrane. 

The simulation domain consists of 200 x 80 lattice units. 
We have tested that when the width is greater than 200 lattice 
units, increase of width dose not affect the simulation results 
if the boundary conditions of both left and right sides are pe- 
riodic. Initially, the fluid was homogeneous and steady. A 
circular membrane placed at the center of the square without 
stretching was discretized into N = 100 segments. The ra- 
dius of the initial circular membrane was 20 lattice units. The 
length in each lattice unit corresponded to 0.14 pia. The re- 
laxation time r was fixed to be 0.75, resulting in 1.69 x 10~^ s 
for each time step. The initial density of the fluid inside and 
outside the close membrane was set to be one lattice Boltz- 
mann unit. The other non-dimensional quantities relevant to 
lattice Boltzmann simulations could be computed correspond- 
ingly. In the simulation, the fluid in the square of 6 x 6 lattice 
units at the center of the system was pumped out with a speed 
of 1.69 X 10^'' g/s, i.e., 1/1000 per time-step, until the prede- 
termined chemical potential drop was reached, the mass inside 
the erythrocyte is then remained constant for the remainder of 
the simulation. 
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Fig. [T] shows the profiles of an erythrocyte with different 
chemical potential drops Afi together with that of a shooting 
method 11911 . As A/i increases, the erythrocyte deforms from 
a circle to an ellipse, and into an biconcave shape. Excellent 
agreement can be found between the two methods. In order 
to further characterize the agreement, we have also computed 
the relative global error <t of the curvature of the membrane 
between the two methods, defined by 



N 



N 



(10) 



i=l 



4=1 



where fc^ and k'^ are the curvatures at segment i, calcu- 
lated from lattice Boltzmann simulations and the shooting 
method lfl9ll . respectively. We found that the relative global 
error was less than 10^^. 




FIG. 2: Static profiles of erythrocytes for bending modulus vary- 
ing according to Eq. i ll lb . (a) For small chemical potential drop 
A/i = 4.24 X 10"'" I/mol, lines with DjOiA, ; correspond to 
n = 2,3,4,5, respectively, (b) and (c) Erythrocyte profiles for larger 
chemical potential drops. Afj. = 0.68, 0.68, 1.87, 1.10, 1.18 and 1.10 
(X lO^^J/mol) for n = 2, 4, 3, 5, 19, 20. 



Due to high values of shear stress in the large arteries or 
in cases of oxidant injury Hf], erythrocyte membranes can be 
pathologically damaged so that the bending rigidity modulus 
becomes non-uniform. To study the effect on the shapes of the 
erythrocytes, we performed numerical simulations with bend- 
ing modulus changing periodically along the membrane 

Kb = Kq [5+{l-S)cos^{m:l')] ,0<l' <1, (11) 

where 5 is a constant, I' is the normalized arc length of the 
membrane, and n is an integer. The simulation results for 
6 = 0.1 are shown in Fig. |2] For small chemical potential 
drop, say Afj, < 4.24 x 10~^ J/mol, the shape of the erythro- 
cyte exhibits the same symmetry of Kb- Remarkably, when 



A/i is large enough, all erythrocytes, with different bending 
modulus, become biconcave shapes, i.e., the biconcave shapes 
can effectively resist the perturbations. We note that the per- 
turbations are quite large as the minimum of Kb is only 10% 
of its original value. Further, erythrocyte profiles for perturba- 
tion wave numbers of the same parity are very similar to each 
other As shown in Fig.|2](b), the difference between the pro- 
files for n = 3 and n = 5, or that between n = 2 and n = 4 
is almost indistinguishable whereas the discrepancy between 
the cases of n = 2 and rt = 3 is clear When the wave number 
n is large enough, the difference for odd and even n becomes 
very small as shown in Fig.|2](c). The chemical potential drop 
needed to collapse an erythrocyte is smaller for an even n than 
that for an odd n, this difference becomes vanishingly small 
as n increases. 



S>30 



1 ' 1 ' r 

(a) 

Green Functions 
.-O- Keller Skalak 2Dy 





•/=598(s ) Slope(Experiment) 
= 0.032 
Slope(LBM_2D) 
= 0.0064 




0.8 

r 



FIG. 3: Erythrocyte deforming and moving in shear flow, (a) Orien- 
tation angle as a function of the swelling ratio F with the shear rate 
7 = 7481s~^, together with the results of Green function and KS 
theory (20ll . (b) Tank tread frequency / plotted against the shear rate 
7. □ and Q correspond to lattice Boltzmann simulation and experi- 
mental results (21I1. respectively. The frequency of lattice Boltzmann 
simulation is multiplied by 27r. Inset: a terminal snapshot of an ery- 
throcyte in a shear flow for 7 — 598 and 449 s^^. 

Now, we performed simulations of an erythrocyte deform- 
ing and moving in a shear flow. According to the KS the- 
ory | |20i| . two parameters relevant to the tank tread- like motion 
of an erythrocyte are the swelling ratio 

r = 47rS'/P^ (12) 



and the viscosity ratio 



^ Pin ^in / Pout ^out ; 



(13) 



where S and P are the area and the girth of a two-dimensional 
erythrocyte, p and v are the density and kinematic viscos- 
ity of the fluid inside or outside the erythrocyte, marked 
by subscripts in and out, respectively. In our simulations 
r « 1 .0 jsl] . At first, we pumped out some fluid inside the ery- 
throcyte until the chemical potential reached a certain value. 
Consequently, the upmost and lowest layer in the simulation 
domain has a rightward and leftward velocity, respectively, 
with a same magnitude to produce a simple shear flow. Mean- 
while, the area of the erythrocyte was kept nearly invariable 
by minus increased mass 

Sp - SqPo 



Ad 



M 



(14) 
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from every fluid node inside the erythrocyte, where Sq and 
S are the areas of the erythrocyte before and after adding a 
shear flow, po ^nd p are the average densities inside the ery- 
throcyte before and after adding a shear flow, M is the number 
of fluid nodes inside the erythrocyte. According to the KS the- 
ory, the orientation angle of an erythrocyte moving in a shear 
flow win vary with F. From Fig. [3] (a), we can see that our 
simulation results agree with the results calculated by Green 
functions lEoll very well when the swelling ratio F > 0.8. In 
fact, when F < 0.8, the shape of the erythrocyte without any 
shear flow will be concave and its terminal shape with a shear 
flow will be protuberant, but different from ellipsoid. In Fig. 
|3](b), the swelling ratio F ~ 0.418, the erythrocytes undergo 
tank tread-like motion, consistent with the famous experimen- 
tal observation by Fischer and Schmidt-Schonbein 112 ill . The 
final shape of the erythrocyte for 7 = 598 s^^ are shown in 
the inset of Fig.[3](b). Numerically we find that the erythro- 
cyte deformed from original biconcave shape into a highly 
flattened shape, which the the boundary integral and singu- 
larity methods can not get despite two-dimensional or three- 
dimensional calculation |0]. We have computed the tank tread 
frequency with respect to the shear rate for F — 0.418. The 
result is displayed in Fig. |3](b) together with the experimen- 
tal one IftIi . It is clear that our two-dimensional simulation 
result has the same linear behavior as the frequency of ex- 
periment. Because the swelling ratio of a cross section of 
three dimensional biconcave erythrocytes may continuously 
increases to one, our two-dimensional simulation is only for 
one certain cross section of three-dimensional, for example 
the center cross section. 

To summaries, we have developed a lattice Boatsman 
model to simulate vesicle deforming and moving in various 
fluid flows for a wide range of Reynolds numbers. The results 
show that the static biconcave shape is rather steady that it can 
effectively resist pathological changes on the membranes. In a 
shear flow, an erythrocyte will deform from biconcave shape 
into a highly flattened one. Considering the intrinsic paral- 
lelity of the lattice boltzmann method, this scheme may have 
good prospects for the simulations of various vesicles deform- 
ing and moving in the vessel or capillary. 
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